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We derive closed analytical expressions for the order parameter $(2;) and for the chemical potential 
pi of a Bose-Einstein Condensate loaded into a harmonically confined, one dimensional optical lattice, 
for sufficiently weak, repulsive or attractive interaction, and not too strong laser intensities. Our 
results are compared with exact numerical calculations in order to map out the range of validity 
of the perturbative analytical approach. We identify parameter values where the optical lattice 
compensates the interaction-induced nonlinearity, such that the condensate ground state coincides 
with a simple, single particle harmonic oscillator wave function. 
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I. INTRODUCTION 

The dynamics of Bose Einstein Condensates (BEC) in 
harmonically confined optical lattices is one of the key 
objects of cold matter research While the fundamen- 
tal equation which describes the mean field solution, the 
Gross-Pitaevski equation (CPE) Q, is readily amenable 
to a numerical solution, explicit analytical expressions for 
the condensate ground state $ and the associated chem- 
ical potential \x are scarce. Variational solutions have 
been proposed in [jj, and quasiclassical approximations 
were implemented, providing several analytical results for 
the CPE Q • Eigenenergies and eigenmodes of the non- 
linear Schrodinger equation in a parabolically confined 
optical lattice can be characterized through asymptotic 
expansions of Mathieu functions, in the tight-binding ap- 
proximation 

We will demonstrate in the present contribution that 
closed analytic results can indeed be derived, and that 
they permit not only qualitative but also quantitative 
insight into the competition between the interaction- 
induced nonlinearity on one hand, and the redistribution 
of the particles over the spatially modulated potential, 
on the other. Comparison with exact numerical results 
will finally allow us to demarcate the range of validity 
of our perturbation approach, with respect to parameter 
values employed in state of the art experiments. 



II. ANALYTICAL EXPRESSIONS 

We consider the one-dimensional nonlinear CPE in a 
stationary optical lattice |2j, 
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where u> is the frequency of the confining harmonic trap 
I™, m the atomic mass, Er — |^;(^f-) 2 the recoil energy 
[7(, s > an adjustable parameter defined by the laser 
intensity, d the laser wavelength, and Aid the atom-atom 
interaction strength. With this one dimensional ansatz, 
we implicitly assume that the dynamics in the transverse 
direction be frozen - e.g., by sufficiently strong confine- 
ment or some other adiabaticity condition [8]. Further- 
more, we imply the number of atoms per site to be large 
enough to justify the mean field limit [Io| . The nor- 
malization condition for $ reads 
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dx\<f>\ 



(2) 



Equation |T]) admits the following integral representa- 
tion for 4>(x) ■[llj]: 
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with the Green function G(x, x') of the linear operator 
Ln = — ^—d 2 ldx 2 + ^muj 2 x 2 , and 

u 2m 1 2 ' 
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the inhomogeneity of the differential equation Lq[$] = f. 

Recall that for Eq. © to hold, we have to impose 
that all involved functions tend to zero as x — > ±00. 
Thus, on the basis of the general theory of Fredholm 
integral equations [H, [H, Qjj], the harmonic oscillator 
functions {<f n } represent a complete basis set which span 
the Hilbert space of the integral equation ([3]) . Using the 
spectral representation of G(x,x'), we can rewrite the 
order parameter $ as 



(4) 



where l a = yj h/muj is the harmonic oscillator wave- 
length. Inserting (H} into ©, the vector coefficient 
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C(u) = {C„(^)}^L is restricted to obey the relation 



[AO) + AC T -T • C-V a P] C = , 



(5) 



with C T the transpose of the vector C. T p i mn and 
Pkm are fourth and second rank tensors, respectively, 
and defined in Appendices A and B. Furthermore, V — 
sEr/Huj, A = Xi£>/l huj, and 



1 _ JL I A 

2 hu } nm 



(6) 



the chemical potential, at second order in A and V : 
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B. Order parameter 



A. Chemical potential 

We tackle the problem under the assumption that the 
atom-atom interaction, the nonlinear term A|$| 2 , and 
the optical lattice potential can be considered as pertur- 
bations with respect to the trap potential mu 2 x 2 /2. In 
this case, the dimensionless chemical potential /i/Huj and 
the vector C in Eq. (O can be expanded in series, 



c = c(°) + c« + c( 2 ) + ... 
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where the quantities and ^ are understood to be of 
the same order in A and V . Expansion of Eq. to sec- 
ond order in A and V leads, after lengthy but straight- 
forward manipulation, the following expression for the 
system's ground state energy: 
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Using the properties of T p i mn and Pfe p given in the Ap- 
pendices A and B we obtain 
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with a = 2ixl j d the ratio of trap to laser wavelength. 
Finally, the series in (|10p can be summed up with the 
help of Eqs. (|CHIC4p in Appendix C, and we obtain for 



Following the same procedure as above, C and thus 
(through ([4])) the normalized order parameter $ can be 
expressed, at first order in A and V , as 

X \ ^ AToo02m — V P()2m ,X 

$ = ^o( r )-.L — Vam(r)- ( 12 ) 
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With the expressions derived in the Appendices A and 
B for the matrix elements T oo2m and Po2mi this can be 
rewritten as 



a = r£^ + V (-i) m+1 J Ay^M! 
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Since, according to Eqs. (|C5IIC8[) of Appendix C, 
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can be condensed into 
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we obtain the following closed expression for the order 
parameter: 
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FIG. 1: (Color online) Variation of the chemical potential 
fi/hu + Vo/2 as a function of the dimensionless self-interaction 
parameter A. Solid, dashed, and dash-dotted lines represent 
the analytical results for d/l = 1, 2 and 3, respectively. Sym- 
bols correspond to the numerical solution of Eq. {l}. The 
strength of the periodic potential was set to V = sEr/Fiw = 
1. Inset: Confining potential U = \mui 2 x 2 — sErcos 2 (^fx) 
for 2 < d/h < 4. 



C. Numerical Solution 

For a numerical solution of Eq. |T]) , we subdivide the 
dimensionless coordinate z = x/l , z G [-L; L], such that 
Zi = (—L/2+(i—l))S, i = 1, . . . , L+l, with a grid size S. 
For 6 <C d/l oy a three-point approximation with uniform 
spacing allows to rewrite the second derivative and the 
differential equation ([TJ as difference equations: 
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where fa = v7o < £ > (•£.;) is the local, normalized order pa- 
rameter which abides by the boundary conditions cf>o = 
4>l+2 = 0. Equation (fT5|) defines a tri-diagonal eigenvalue 
problem with the normalization condition l^l 2 ^ = 1 
or, equivalently, ^ i \4>i\ 2 — 1, where fa = Since 
the resulting matrix (| 1 5|) involves the non-linear term 

| fa | 2 , an iterative procedure is needed for a self-consistent 
solution. To avoid numerical instabilities [lil . [l6| , we in- 



FIG. 2: Same as in Fig. [TJ for V a = 1 (solid line), 30 (dashed 
line), 60 (dash-dotted line), and 200 (dotted line). Symbols 
represent the analytical results. The laser wavelength was set 
to d/l =0.4. 



serted the iterative wave function 
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in (|15p . at each new iteration fc, with e£ [0; 1], and start- 



ing condition F t 



(0) 



This procedure was repeated 



until 



< 8+ and|(^( fc ) - ^ k -V)/hLu\ < 5^ 

where 6^ and <5 M are the desired accuracies for the order 
parameter and the chemical potential, respectively. Our 
iterative solutions then converge to a relative uncertainty 
of 10 -4 for fi/huj and the normalized order parameter. 

III. RESULTS 

Figure [Tj displays the variation of the chemical po- 
tential, n/huj + 0.5V o , as a function of the dimension- 
less atom-atom interaction A. The analytical solution 
dTTJ|) is represented by solid, dashed, and dash-dotted 
lines, while symbols represent the numerical solutions 
of Eq. JTJ. The inset illustrates the confining potential 
U = rnuj 2 x 2 /2 — sErcos 2 (2irx/d), as a function ofx/l , 
for V = 1 and 2 < d/l < 4. In our calculation of fi, we 
fixed V — 1 and checked the variation of the perturbative 
solution as compared to the numerical one, for several 
values d/l a = 1,3 and 4 of the normalized wavelength. 
For A < the relative error between the numerical so- 
lution and Eq. (fTTj) increases as d/l increases, while for 
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(a) A=0.0 




(b) A= -0.5 




(d) A=-1.5. 
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FIG. 3: Normalized order parameter ylo<&(x/lo) for the di- 
mensionless self-interaction values: (a) A =0, (b) -0.5, (c) -1.0, 
and (d) -1.5. V — 1 and d/l =1. Solid line: Closed solution 
according to Eq. (|14p . Open circles: Numerical result. 



A > the worst comparison is obtained at d/l a = 1. In 
the case of repulsive interaction, numerical and pertur- 
bative calculation exhibit small variations with d/l . In 
general, Eq. (fTTjl describes the values of the chemical po- 
tential better for repulsive than for attractive interaction. 
This can be understood with the following argument: In 
[Til [Llj it was shown shown that, at V — 0, pertur- 
bation theory matches the numerical solution very well 
in the interval — 2 < A < 2. When the optical lattice 
is switched on (V ^ 0), we can argue that an effective 
renormalization of the non-linear parameter takes place. 
Since Eq. (fTTj) includes a negative quadratic term in V a , 
we expect that the effective range of A where perturba- 
tion theory is applicable will shift to A > —2, i.e., more 
repulsive or less attractive atom-atom interactions. Nat- 
urally, the behavior of the chemical potential with d/l a 
limits this simple description. 

In Fig. [51 we plot the dependence of /i(A)/?kj + 0.5V 
on A, for several values of the strength of the periodic 
potential V . In the calculation we took d/l Q = 0.4. The 
figure shows perfect agreement of Eq. (fTT|) with the nu- 
merical result for V ranging between 1 and 60, while, at 
V = 200, we observe large discrepancies. 

Figure [3] displays the normalized order parameter, 
&(x/lo)Vki> f° r f° ur values of A = 0,-0.5,-1.0, and 
— 1.5, and fixed V — 10, d/l Q = 1. In general, we 



FIG. 4: Order parameter ®(x/lo) for d/l a = 1. (a) A =2, 
Va = sEr/Klo = 0.5, 1, 5, and 10. (b) V a = 10, A =-2, 0, 
2. Symbols represent the numerical solution. For the sake of 
comparison, U (see caption Fig. [TJ and the harmonic oscillator 
potential are shown by solid and dashed lines, respectively. 
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obtain a remarkably good agreement between the ana- 
lytical prediction (fT4")l and the numerical result, for all 
considered (attractive) values of the dimensionless inter- 
action parameter A, and witness the interplay between 
the optical lattice and harmonic oscillator potential: The 
harmonic oscillator ground state tp is modulated by an 
oscillatory behavior, induced by the optical lattice po- 
tential. We further note that the wave function gets 
more localized, and the maximum of §(x/Iq) increases, as 
the non-linear potential becomes more attractive. Panel 
(a) of Fig.[l]shows the variation of $(x/io)\/^o with the 
strength of the periodic potential V , for the repulsive 
case A = 2. As V increases from 0.5 to 10, the lattice- 
induced density oscillations manifest with increasing am- 
plitude. A comparison between attractive and repulsive 
interactions, for V — 10, is given by panel (b) of Fig. [4] 
Since d/l = 1 in Figs.[3]and|ll the wave function $(x/io) 
exhibits oscillations at x/Iq w n/2 (n — ±1,±2, ...), 
which are quenched because of the asymptotic decay oc 
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FIG. 5: (Color online) Quenching of the non-linear repulsive 
interaction. Left panel: Region Ti is defined by the condition 
(fT7)l . The curve V is defined by Eq. (JT7J, at x = 0. Right 
panel: Order parameter calculated for those parameter values 
identified by A, B, and C on the left. Solid lines represent $ 
as derived from Eq. (|14| , while open circles indicate the exact 
numerical solution. Dashed lines correspond to the harmonic 
oscillator wave function. 
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FIG. 6: (a) Range of validity of Eq. JTT| in the V -d/l - 
plane. The functions Linear, Mixed and Quadratic (see text) 
are represented by dashed, dotted and solid lines, respectively. 
z c = [(Vo)c, (d/lo) c ] is defined by equality between Linear 
and Quadratic, (b) Zoom into the region of large values of 
V - The circle locates theparameter values employed in the 
experiments reported in (171 ]. 



exp (— x 2 /2ll) of the oscillator wave functions. These os- 
cillations are illustrated in Fig.[2Jb), where the solid line 
represents the potential U (see inset of Fig. [T]), along with 
the wavefunction $ for the same values of V and d/l Q . 
The positions of the observed maxima (minima) in the 
order parameter correspond to the minima (maxima) of 
the resulting effective potential U. 

The extrema of $ exhibit different behavior depend- 
ing on the non-linear interaction A. As the condensate 
becomes effectively less confined for increasing repulsive 
interaction, the amplitudes of the maxima and minima 
of $ are enhanced. From the physical point of view it is 
clear that a strongly attractive interaction (in the figure: 
A = —2) should flatten the lattice-induced oscillations 
of the condensate. In panel (b) we compare the evolu- 
tion of the order parameter from attractive to repulsive 
interaction. We note that, as A increases, the conden- 
sate spreads, i.e., the wave function delocalizes, while its 
maximum decreases. The analytical solution (|14[) is less 



accurate at x ~ and A = — 2 (with an error less than 
5%) , in agreement with the results of Fig. [1] and the an- 
ticipated range of validity of perturbation theory. Nev- 
ertheless, in this particular case the observed agreement 
with the exact solution is quite remarkable at x ^ 0, and 
shows that the closed analytical expressions derived here 
allow a fair representation of the BEC wave function in 
an optical lattice. Almost no differences can be observed 
on the scale of the figure between numerical and analyt- 
ical solutions for A = 2. 



A. Quenching of the non-linear interaction 

Since the repulsive atom-atom interaction and the opti- 
cal lattice potential have opposite signs, the set of param- 
eters (Ai£>, u>, d, sEr) can be chosen such as to minimize 
the effect of the non-linearity. To study this quenching 
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effect, in the vaiidity range of our perturbation approach, 
we consider the order parameter given by Eq. HU). In 
the inset of the left panel of Fig. 03 we monitor the de- 
pendence of G{x/l ) and F(x/l ,a) on x/l a . While the 
function G is universal, F shows a strong dependence on 
the parameter djl . Moreover, for d/l = 4 the contribu- 
tion of F to the order parameter, induced by the optical 
lattice, resembles the function G, though with opposite 
sign. If we impose the condition 

AG(x/l o ) + V o F(x/l o ,a)^0, (17) 

then $ can be solely described by the harmonic oscilla- 
tor wave function ip . The range of values of A/V a , d/l , 
and x/l Q that satisfy (fT7|) is represented by the hatched 
region H in Fig. [5] At x = 0, Eq. (fl"7|) defines the curve 
T> shown in the figure. We have noted that at x = the 
functions G and F reach their minimum and maximum 
values, respectively. In the right panel of the figure we 
compare the harmonic oscillator wave function ip with 
the order parameter as predicted by (|14p. We chose three 
different points A, B, and C in the parameter space of 
the left panel for this comparison: Indeed, A does not 
fulfill the condition (jTTJ) , and, consequently, <f> does not 
match tp . The case C belongs to region H, but not to 
the curve T>. Hence, we observe a small discrepancy be- 
tween both functions, mainly at x — 0. Case B, which 
sits right on top of the curve T>, yields a perfect match 
between harmonic oscillator wave function and order pa- 
rameter. Exact numerical solution of the GPE (Q]) with 
the parameters defined by B corroborates this result, as 
indicated by the open circles in the figure. 



B. Validity of the perturbation approach 

For a vanishing lattice, V — 0, our closed analytical so- 
lution for the chemical potential and its comparison with 
the numerical solution provides universal criteria for its 
range of validity: in the interval |A| < 2, we derive an 
accuracy better than 97% [H for Eq. (jnj. For non- 
trivial values of the lattice strength, two more indepen- 
dent parameters, V and d/l Q come into game. A neces- 
sary condition for the validity of the perturbative results 
requires that the functions Linear = V exp (— a 2 ) /2, 

Mixed = AV exp (-a 2 ) {J5i(f ) - C - lnf } /v^tt, and 
Quadratic = V 2 exp (-2a 2 ) {Chi(2a) - C - In 2a} /4, 
that contribute to Eq. (fTTj) . simultaneously be smaller 
than one. This criterion allows to confine the range of 
validity of our analytical results in a 2D map spanned by 
V and d/l , as displayed in Figs. [Ha) and (b). There 
the Mixed contribution to (jTTJ) is plotted for three dif- 
ferent values of |A| = 0.5, 1, and 2. Fig. thus de- 
fines different V vs. d/l Q regions where Eq. (JTTJ) can be 
trusted: In region I, the quadratic term in V is less im- 
portant than the linear one, while, in region II, we need to 
make one more distinction, depending on the d/l Q values. 
For example, if d/l a < (d/l ) c , we have that Quadratic 



represents the main contribution to /i, and vice versa 
if d/l > (d/l ) c , where Linear > Quadratic. A similar 
analysis applies in regions III and IV, with respect to the 
relative contributions of Linear and Mixed. 

Finally, let us note that experimental data as reported 
in [17j on S7 Rb condensates were obtained for V = 275 
and d/l = 0.19. This working point is indicated by a 
solid circle in Fig. [^b), and shows that our perturba- 
tive method is applicable in the experimental parameter 
range. 



IV. CONCLUSIONS 

We provided closed analytical expressions for the 
chemical potential and the particle density of the GPE 
ground state in an optical lattice. These solutions 
were obtained by a perturbative expansion in the lattice 
strength and in the (attractive or repulsive) atom-atom 
interaction, with a well-controlled range of validity in the 
associated parameter space. 

Interestingly, the interaction-induced non-linearity 
may be quenched by the presence of the lattice. Under 
certain conditions (see left panel of Fig. [5]), we predict 
that the particle density of the repulsively interacting sys- 
tem is given by $(x,t) = exp(i/z[A, V , d/lo]t/h)tp (x/lo). 

The solutions derived here can serve as a useful tool to 
study a weakly interacting BEC in a not too deep ID op- 
tical lattice. On this basis, it is possible to develop a uni- 
fied and comprehensive picture of Bogoliubov equations, 
the time-dependent GPE and of collective excitations (a 
subject under investigation). The model here developed 
can be generalized to two and three dimensions, and also 
for two component BECs. 
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APPENDIX A: TENSOR T 

Using the harmonic oscillator wave function 

vnix) = 0^4^) 7 exp Hn (£) ' 

. (Al) 

with H n (z) the Hermite polynomials, the matrix ele- 
ments Tpi mn are defined by 
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It is possible to perform the above integral, and we get 

Mm 
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T(M -l + l)T(M -n + \) 
Y(M-n-l + \) 



3F2 



/ -m, -p, -M + n + l + \ ■ \ 
\ -M + Z + ±, -M + n+\; ) ' 



(A3) 

with T(z) the gamma function, 3F2 ' a ' a ' 1 the 

\ Pi, P2, / 

generalized hypergeometric function [l9| , and 2M = p + 
Z + m + n. Consequently, we arrive at the following useful 
relations: i) T p i mn = 0, if n + m + l +p is an odd number, 
ii), for to = 0, 



TplnO — 



2 s ' 1 T(s - l)T(s - p)T(s - n) 



y/2 n + l +H\p\n\ 
and, iii), for p = I = n = 0, 
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APPENDIX B: TENSOR P 

The contribution of the optical lattice is represented 
by the two dimensional matrix elements P kp given by 



Pk. 



1 



hp 



y/n2 k +Pk\p\ J- 



,2-kI, 



cos (— — z) X 



H k (z)H p (z) exp(-z 2 )] dz. 
For p = k + 2m this integral is equal to pj 



(Bl) 



1 



hp 



2 m+1 V (k + 2m)! 



m.O 



(-l) m 6 2m exp(-fo 2 /4)i 



2)>! 



, (B2) 



where b = Anl /d and Li (x) are the Laguerre polynomi- 
als. The symmetry of the Hermite polynomials imposes 



that P kp = if p = k + 2m + 1. Using (|B2[) the following 
relations hold: i) For k = p = 0, 



1 1 

p oa = 2 + 2 exp 



and, ii), for k = and p = 2to, 
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APPENDIX C: SERIES 



We can sum up the series 



Fl 



3_ (2m- 1)! 
~~ 2 4 '»(to!) 2 

71=1 V ' 
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noting that 
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(2m- 1) 
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1 2x 2 + l 



2 V 2(x 2 + 1) 



Hence, Fl = 0.033106. 

Furthermore, we have the series [l| 
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■a 4m = Chi(2a)-C-\n2a, (C4) 



^ 2to(2to) 

where Ei(x) is the exponential integral, Chi[x) is the 
cosine hyperbolic integral, and C is Euler's constant. 

For the summation of the first part of the series in (| 13[) , 
we note that 



F2(z,a) 



~ (~2) m a 2m 

^1 V( 2m ) ! 

cxp(-z 2 /2) 



2\/W 72 



(exp a 2 cos 2az — l) . (C5) 



The series 
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is related to the function F2(z, a) through the differential Now the series 
equation 



dGl(z,a) _ F2{z,a) 
da a 

The solution of l|C7[) is given by 



(C7) 
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m2 m m\ 



(CIO) 

is related to the function F3(z, c) through the equation 



Gl(z, a) 



1 / z 

exp - — 



2v/^rV2 

- [exp(y 2 )coa(2yz)-l]dy. (C8) 
o V 



dg 2(3 



dfi (1-/3 2 ) 



F3. 



(Cll) 



To perform the sum of the first series in CEO), we use that with P = 1 / c - The differential equation JU11J) admits the 

solution 



F3(z,c) 



1 



cexp(-c z + — ) - exp - — 



2- 2^! I 1 "?' lP2m[z) - 

m—1 x 



g (,c )= 2exp i3) r exp( " g(1 :/ 2)) "\. 



(C9) 



(C12) 
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